* Figure VII panels e)-f)
* 29 Oct 2021

clear all
set more off

grstyle init
grstyle set plain, nogrid horizontal
grstyle color background white

*Choose figure width:
local fig_w= 5
local res = `fig_w' * 600

import excel "$output_s/w5_test_of_model.xlsx", sheet("Sheet1") firstrow clear 

* Bar graph of cases - model vs prediction
gen occup_opt = ceil(case_opt_w5/2) //Quick way to separate cases
gen occup = ceil(case_real_w5/2)

label define occupation 1 "Mixing" 2 "Livestock Rearing Only" 3 "Wage Work Only" 4 "No Work"
label values occup* occupation

drop if occup == 4

preserve
	keep hhid occup*
	rename occup_opt occup2
	rename occup occup1
	reshape long occup, i(hhid) j(tag)
	label define tags 1 "Actual" 2 "Predicted"
	label values tag tags
	catplot tag occup, percent(tag) var1opts(label(labsize(small))) var2opts(label(labsize(small))) recast(bar) asyvars bargap(20) bar(1, color(pink*0.4)) bar(2, color(green*0.4)) ytitle("Percentage") name(bars, replace) yla(0(10)50, gmax)
	graph save bars "$output/Bars_cases_excl_no_work_w5.gph", replace
	graph export "$output/Bars_cases_excl_no_work_w5.png", as(png) replace
restore


set level 95 //Setting confidence interval


* Figure VII panel e)

	*Colored version:
	twoway (lpolyci l_opt_w5 k_w5 if k_w5<16000, lcolor(green) fcolor(green%30) alcolor(%30) bw(800)) (lpolyci l_w5 k_w5 if k_w5<16000, lcolor(pink) fcolor(pink%30) alcolor(%30) bw(800)), xtitle("Capital at Wave 5 (BDT)", margin(t+2)) ytitle("Hours in Livestock Rearing") yla(0(200)800) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_l, replace)
	graph save Lpoly_l "$output/Lpoly_l_kw5.gph", replace
	graph export "$output/Lpoly_l_kw5.png", as(png) replace
	
	*TIFF Greyscale version:
	twoway (lpolyci l_opt_w5 k_w5 if k_w5<16000, lcolor(gs10) lp(dash) fcolor(gs10%30) alcolor(%30) bw(800)) (lpolyci l_w5 k_w5 if k_w5<16000, lcolor(gs3) fcolor(gs3%30) alcolor(%30) bw(800)), xtitle("Capital at Wave 5 (BDT)", margin(t+2)) ytitle("Hours in Livestock Rearing") yla(0(200)800) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_l, replace) xsize(`fig_w')
	graph export "$output/TIFF/Lpoly_l_kw5.tif", as(tif) replace width(`res')

* Figure VII panel f)

	*Colored version:
	twoway (lpolyci h_opt_w5 k_w5 if k_w5<16000, lcolor(green) fcolor(green%30) alcolor(%30) bw(800)) (lpolyci h_w5 k_w5 if k_w5<16000, lcolor(pink) fcolor(pink%30) alcolor(%30) bw(800)), xtitle("Capital at Wave 5 (BDT)", margin(t+2)) ytitle("Hours in Wage Work") yla(0(200)1400) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_h, replace)
	graph save Lpoly_h "$output/Lpoly_h_kw5.gph", replace
	graph export "$output/Lpoly_h_kw5.png", as(png) replace

	*TIFF Greyscale version:
	twoway (lpolyci h_opt_w5 k_w5 if k_w5<16000, lcolor(gs10) lp(dash) fcolor(gs10%30) alcolor(%30) bw(800)) (lpolyci h_w5 k_w5 if k_w5<16000, lcolor(gs3) fcolor(gs3%30) alcolor(%30)  bw(800)), xtitle("Capital at Wave 5 (BDT)", margin(t+2)) ytitle("Hours in Wage Work") yla(0(200)1400) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_h, replace) xsize(`fig_w')
	graph export "$output/TIFF/Lpoly_h_kw5.tif", as(tif) replace width(`res')
	